
typedef double x_t ;
typedef double real_t ;

void vl_imwbackward_d
(x_t const* src_pt, int M, int N,
 real_t * x_pt, real_t * y_pt,
 
 x_t *dst_pt, int K)
{
                      

  X_pt = getPr(X); 
  Y_pt = getPr(Y) ;
  I_pt = getPr(I) ;
  iwXp_pt = getPr(iwXp) ;
  iwYp_pt = getPr(iwYp) ;

  /* Allocate the result. */
  out[wI] = mxCreateDoubleMatrix(Mp, Np, mxREAL) ;
  wI_pt = mxGetPr(out[wI]) ;

  /* -----------------------------------------------------------------
   *                                                        Do the job
   * -------------------------------------------------------------- */
  Xmin = X_pt[0] ;
  Xmax = X_pt[N-1] ;
  Ymin = Y_pt[0] ;
  Ymax = Y_pt[M-1] ;

  for(jp = 0 ; jp < Np ; ++jp) {
    for(ip = 0 ; ip < Mp ; ++ip) {    
      /* Search for the four neighbors of the backprojected point. */     
      double x = *iwXp_pt++ ;
      double y = *iwYp_pt++ ;
      double z = NaN ;

      /* This messy code allows the identity transformation
       * to be processed as expected. */
      if(x >= Xmin && x <= Xmax &&
         y >= Ymin && y <= Ymax) {
        int j = findNeighbor(x, X_pt, N) ;
        int i = findNeighbor(y, Y_pt, M) ;
        double* pt  = I_pt + j*M + i ;
        
        /* Weights. */
        double x0 = X_pt[j] ;
        double x1 = X_pt[j+1] ;
        double y0 = Y_pt[i] ;
        double y1 = Y_pt[i+1] ;
        double wx = (x-x0)/(x1-x0) ;
        double wy = (y-y0)/(y1-y0) ;

        /* Load all possible neighbors. */
        double z00 = 0.0 ;
        double z10 = 0.0 ;
        double z01 = 0.0 ;
        double z11 = 0.0 ;

        /* 2006: possible bug? should pt++ in any case? */
        if(j > -1) {
          if(i > -1 ) z00 = *pt ;
          pt++ ;
          if(i < M-1) z10 = *pt ;
        }

        pt += M - 1;

        if(j < N-1) {
          if(i > -1 ) z01 = *pt ; 
          pt++ ;
          if(i < M-1) z11 = *pt ;
        }       

        /* Bilinear interpolation. */
        z = 
          (1-wy)*( (1-wx) * z00 + wx * z01) +
              wy*( (1-wx) * z10 + wx * z11) ; 
      }
      
      *(wI_pt + jp*Mp + ip) = z ;
    }
  }
}
